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ABSTRACT 

The BLAST programs are widely used tools for 
searching protein and DNA databases for sequence 
similarities. For protein comparisons, a variety of 
definitional algorithmic and statistical refinements 
described here pemiits the execution time of the 
fnhi •'^^°?u®.'"® *° decreased substantially while 
enhancing their sensitivity to weak simliarltles. A new 
criterion for tnggering the extension of word hits 
combined with a new heuristic for generating gapped 

a KorolS; ^l^^^J S^fP^'' ^"^2^ P'°3ram that mns 
fn ^PP^fx^ately three times the speed of the original, 
in addition, a method is introduced for automatically 

5n!?rtK!?i.^?"?^''y s'9"!«cant alignments pro- 
duced by BLAST into a position-specific score matrix, 

at A^-^^ Position-Specific Iterated BLAST (PSI- 
LiS^n<frT♦^'■^^ '""^ approximately the same 
l^t^' ^® S^PP®^ BLAST, but in many 

cases is much more sensitive to weak but biologically 
relevant sequence similarities. PSI-BLAST is used to 

sloTsu^rS^ir '"'"''''^ °' 

INTRODUCTION 

Variations of the BLAST algorithm (1) have been incorporated 
mto several popular programs for searching protein and DNA 
oatabases for sequence similarities. BLAST programs have been 
witten to compare protein or DNA queries with protein or DNA 
dajbases m any combination, with DNA sequences often 
undergoing conceptual translation before any comparison is 
performed. We will use the blastp program, which compares 
protem queries to protein databases, as a prototype for BLAST, 
auhough the Ideas presented extend immediately to other 
versions that involve the translation of a DNA query or database. 
Sa ^It^^ refinements described are applicable as weU to 
^iNA-DNA comparison, but have yet to be implemented 



BLAST IS a heuristic that attempts to optimize a specific 
smulanty measure. It permits a tradeoff between speed and 
sensmvity, witii the setting of a 'threshold' parameter, T. Ahiaher 
value of T yields greater speed, but also an increased probabiUty 
of mjssmg weak similarities. ITie BLAST program requires time 
prqwnional to the product of the lengths of the query sequence 
and the database searched. Since the rate of change in database 
sizes cunrently exceeds that of processor speeds, computers 
niniung BLAST are subjected to increasing load. However the 
conjuncuon of several new algorithmic ideas allow a new version 
of BLAST to achieve improved sensitivity at substantially 
augnented speed. -Qus paper describes three major fefinemena 
to oLAST. 

(i) For increased speed, the criterion for extending word pain 
has been modified. The original BLAST program seeks short 
word pairs whose aligned score is at least T. Each such 'hit' is then 
extended, to test whether it is contained within a high-scoring 
ahgnment For the default T value, this extension step consume! 
inost of the processing time. The new 'two-hit' method requires 
flie existence of two non-overiapping word pain on the same 
diagonal, and within a distance A of one another, before an 
extension is invoked. To achieve comparable sensitivity, the 
threshold parameter T must be lowered, yielding more hits than 
previously. However, because only a small fraction of ^feese hits 
decreSes"^^ *e average amount of computation required 

(ii) TTie ability to generate gapped alignments has been added 
me ongmal BLAST program often finds several alignments 
mvolvmg a smgle database sequence which, when considered 
tofether, are statistically significant Overiooking any one of 
these alignments can compromise the combined result By 
mtroducmg an algorithm for generating gapped alignments, it 
becomes necessary to find only one rather than all the ungapped 
alignments subsumed in a significant result This allows the T 
parameter to be raised, increasing the speed of the initial database 
scan. The new gapped alignment algorithm uses dynamic 
programmmg to extend a central pair of aligned residues in both 
directions. For speed, earher heuristic methods (2,3) confmed the 
alignments produced to a predefmed strip of the dynamic 
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programming path graph (4). Our approach considers only 
alignments that drop in score no more ^an below the best 
score yet seen. The algorithm is able thereby to adapt the region 
of the path graph it explores to the data. 

(iii) BLAST searches may be iterated, with a position-specific 
score matrix generated from significant alignments found in 
round / used for round / + 1. Motif or profile search methods 
frequently are much more sensitive than pairwise comparison 
meAods at detecting distant relationships. However, creating a 
set of motifs or a profile that describes a protein family, and 
searching a database with them, typically has involved running 
several different programs, with substantial user intervention at 
various stages. The BLAST algorithm is easily generalized to use 
an arbitrary position-specific score matrix in place of a query 
sequence and associated substitution matrix. Accordingly, we 
have automated the procedure of generating such a matrix from 
the output produced by a BLAST search, and adapted the BLAST 
algorithm to take this matrix as input The resulting Position- 
Specific Iterated BLAST, or PSI-BLAST, program may not be as 
sensitive as the best available motif search programs, but its speed 
and ease of operation can bring the power of these methods into 
more common use. 

After describing these refinements to BLAST in greater detail, 
we consider several biological examples for which the sensitivity 
and speed of die program are greatly enhanced 

STATISTICAL PRELIMINARIES 

To analyze the BLAST algorithm and its refinements, we need 
first to review the statistics of high-scoring local alignments. 
BLAST employs a substitution matrix, which specifies a score 
for aligning each pair of amino acids i and / Given two sequences 
to compare, the original BLAST program seeks equal-length 
segments within each that, when aligned to one another without 
gaps, have maximal aggregate score. Not only the single best 
segment pair may be found, but also other locally optimal pairs 
(3,5-7), whose scores cannot be improved by extension or 
trimming. Such locally optimal alignments are caJled *high-scor- 
ing segment pairs' or HSPs. 

For the sake of the statistical theory, we assume a simple protein 
model in which the amino acids occur randomly at all positions 
with background probabilities We require that the expected 

score for two random amino acids ^^j^j^ij be negative. Given 

u 

the Pi and ^y, the basic theoiy (8,9) yields two calculable 
parameters, X and K, which can be used to convert nominal HSP 
scores to normalized scores, thereby rendering all scoring 
systems direcdy comparable from a statistical perspective. The 
normalized score 5' for an HSP is given by the equation: 



In this paper, a nominal score is given without units, while a score 
normalized by equation 1 is said to be expressed in bits (10,11). 
When two random protein sequences of sufficient lengths m and 
n are compared, the number E of distinct HSPs with normalized 
score at least S' expected to occur by chance is well approximated 
by: 

E = N/2^' 2 



where W = mn is the search space size (8-10). If a protein is 
compared to a whole database rather than a single sequence, n is 
the database length in residues. Equation 2 may be inverted to 
yield S' = log2(WE), the normalized score required to achieve a 
particular £- value. In a typical current database search, a protein 
of length 250 might be compared to a protein database of 50 000 OOU 
total residues. To achieve a marginally significant £-value of 
0.05, a normalized score of -38 bits is necessary. 

While the theory just outiined has not been proved for gapped 
local alignments and their associated scores, computational 
experiments strongly suggest that it remains valid (3 ,12-15). The 
statistical parameters X and K, however, are no longer supplied by 
theory but must be estimated using comparisons of either 
simulated or real but unrelated sequences. To distinguish below 
whether a given set of parameters X 2nd K refer to gapped or 
ungapped alignments, we use the subscripts g and u respectively. 

When gaps are not allowed, a further important theorem states 
that within HSPs the aligned pair of letters (ij) tends to occur with 
the *target frequency': 

^.=P^/-5u 3 

The ^ij of equation 3 sum to 1; indeed, 7^ is calculated to be the 
unique positive number for which this is the case (8 ,9). The scores 
5ij are optimal for detecting alignments with these particular target 
frequencies (8,10), and by inverting equation 3 to sij = 
Au, scores may be chosen, with arbitrary scale, that 
correspond to any desired set of qiy The popular PAM ( 1 6, 1 7) and 
BLOSUM ( 1 8) substitution matrices are constructed with expHcit 
use of this log-odds formula. No corresponding result has been 
. established for gapped alignment scoring systems. However, if 
the gap costs used are sufficientiy large, it is expected that the 
target frequencies observed in high-scoring local alignments of 
random sequences will not differ greaUy from those for the 
no-gap case. 

REFINEMENT OF THE BASIC ALGORITHM: THE 
TWO-HIT METHOD 

The central idea of the BLAST algorithm is that a statistically 
significant alignment is likely to contain a high-scoring pair of 
aligned words. BLAST first scans the database for words 
(typically of length three for proteins) that score at least T when 
aligned with some word within the query sequence. Any aligned 
word pair satisfying this condition is called a hit. Th^econd ';tep 
of the algorithm checks whether each hit lies within an alignment 
with score sufficient to be reported. This is done by extending a 
hit in both directions, until the running alignment's score has 
dropped more thanX below the maximum score yet attained. This 
^ extension step is computationally quite costiy; with the T and X 
parameters, necessary to attain reasonable sensitivity to weak 
alignments, the extension step typically accounts for >90% of 
BLAST'S execution time. It is therefore desirable to reduce the 
number of extensions performed 

Our refined algorithm is based upon the observation tliat an 
HSP of interest is much longer than a single word pair, and may 
therefore entail multiple hits on the same diagonal and within a 
relatively short distance of one another. (The diagonal of a hit 
involving words starting at positions {X},X2) of the database and 
query sequences may be defined as xi -xi- The distance between 
two hits on the same diagonal is the difference between their fiJ^^ 
coordinates.) This signature may be used to locate HSPs more 
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Figure L The sensitivity of the two-hit and one-hit heuristics as a function of 
HSP score. Using the BLOSUM-62 amino acid substimtion matrix (18). and the 
target fiequcpcies q\j implied by equation 3 and the background amino acid 
frequencies P\ of Robinson and Robinson (20), 100 000 model HSPs were 
generated for each of the nominal scores 37-92. conesponding to nonnalized 
scores 1 9.9-45. 1 bits. It was determined by inspection whether each HSP failed 
to contain two non-overlapping lengih-3 word pairs with nominal score at le&st 
11 , and within a distance 40 of oiie another, or a single length-3 woid pair whh 
nominal score at least 13. The corresponding probabilities of missing an HSP 
using the iwo-hit heuristic with r= 11, and the one-hit heuristic witii 7= 13. 
are plotted as a function of normalized HSP score. The two-hit method is more* 
sensitive for HSPs with score at least 33 bits. 



efficiently. Specifically, we choose a window length A, and 
invoke an extension only when two non-overlapping hits are 
found within distance A of one another on the same diagonal. Any 
hit that overlaps the most recent one is ignored. Efficient 
execution requires an array to record, for each diagonal, the first 
coordinate of the most recent hit found Since database sequences 
are scanned sequentially, this coordinate always increases for 
successive hits. The idea of seeking multiple hits on the same 
diagonal was first used in the context of biological database 
searches by Wilbur and Lipman (19). 

Because we require two hits rather than one to invoke an 
extension, the threshold parameter T must be lowered to retain 
comparable sensitivity. The effect is that many more single hits 
are found, but only a small fraction have an associated second hit 
on the same diagonal that triggers an extension: The great 
majority of hits may be dismissed after the minor calculation of 
looking up, for the appropriate diagonal, the coordinate of the 
most recent hit, checking whether it is within distance A of the 
current hit's coordinate, and finaOy replacing the old wifli the new 
coordinate. Empirically, the computation saved by requiring 
fewer extensions more than offsets the extra computation 
required to process the larger number of hits. 

To smdy the relative abilities of the one-hit and two-hit methods 
|o detect HSPs of varying score, we model proteins using the 
background amino acid frequencies of Robinson and Robinson 
v20), and use the BLOSUM-62 substitution matrix (18) for 
sequence comparison. Given these Pi and s\y the statistical 
parameters for ungapped local alignments are calculated to be 
^ = 0.3176 and = 0.134. Using equation 3 above, we may 
calculate the q\y for which the scoring system is optimized, and 
employ these target frequencies to generate model HSPs. Finally, 
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Figure 2. TTie BLAST comparison of broad bean leghemogiobin I (87) 
(SWISS-PROT accession no. P02232) and horse ^-globin (88) (SWISS-PROT 
accession no. P02062), Tlie 15 hits with score at least 13 are indicated by plus 
signs. An additional 22 non-oveiiapping hits with score at least 11 are indicated 
by dots. Of these 37 hits, only flie two indicated pairs are on the same diagonal 
and wiflun distance 40 of one another. Thus the two-hii heuristic with 7=11 
triggers two extensions, in place of the 15 extensions invoked by the one-hit 
heuristic with T- 13. Because this is just one example, the relative numbers of 
hits and extensions at the various settings of T correspond only roughly to the 
ratios found in a full database search. An ungapped extension of the leftward 
of the two hit pairs yields an HSP with nominal score 45, or 23.6 bits, calculated 
using ^ and K^^. 



we evaluate the sensitivity of the one-hit and two-hit BLAST 
heuristics using these HSPs. 

The one-hit method will detect an HSP if it somewhere contains 
a length- W word of score at least T. For 1^= 3 and r= 13, Figure 
1 shows the empirically estimated probability that an HSP is 
missed by this method, as a function of its normalized score. The 
two-hit method will detect an HSP if it contains two non- 
overlapping length-W' words of score at least T, with starting 
positions that differ by no more than^i residues. For W^=3, r= 1 1 
and i4 = 40, Figure 1 shows the estimated probability that an HSP 
is missed by this method, as a function of its normalized score. For 
HSPs with score at least 33 bits, the two-hit heuristic is more 
sensitive. 

To analyze the relative speeds of the one-hit and two-hit 
methods, using the parameters studied above, we note that the 
two-hit method generates on average -3.2 times as many hits, but 
only jO. 14 times as many hit extensions (Fig. 2). Because it takes 
approximately one ninth as long to decide whether a hit need be 
extended as actually to extend it, the hit-processing component of 
the two-hit method is approximately twice as rapid as the same 
component of the one-hit method. 

TRIGGERING THE GENERATION OF GAPPED 
ALIGNMENTS 

Figure 1 shows that even when using the original one-hit method 
with threshold parameter 7=13, there is generally no greater than 
a 4% chance of nndssing an HSP with score >38 bits. While this 
would appear sufficient for most purposes, the one-hit default T 
parameter has typically been set as low as 11, yielding an 
execution time nearly three times that for T = 13. Why pay this 
price for what appears at best marginal gains in sensitivity? The 
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Leghenoglobin 43 FSFLKDSAGWDSPKLGAHASKVFGMVRDSAVQLRATGCW--LD6KDGS 90 

P t + V+ ♦PKt- AH -HCV L + 6E V LD 

Beta globin 45 PGDLSKPGAVMGHPKVKAHGKKV tHSFGEGVHHLDHLKSTFAALSE 90 

Leghenoglobin 91 ZKXQKGVLDP-HFWVKEALLKTIKEASGDKW5EELSAAWEVAYD6IATAI 140 

+H K +DP +F ++ L+ + C +f EL A-^-m- G*A A+ 
Beta globin 91 LHCOKLRVDPENFELLCimAnrVLARKFGKDFTPELQASrQKVVAGVAKAL 141 



Figure 3. A gapped extension generated by BLAST for die comparison of broad bean leghemoglobin I (87) and horse ^globin (88). (a) The region of the path graph 
explored when seeded by the alignment of alanine residues at respective positions 60 and 62. This seed derives from the HSP generated by the leftward of the two 
ungapped extensions illustrated in Figure 2. The Xg dropoff parameter is the nominal score 40, used in conjunction with BLOSUM-62 substitution scores and a cost 
of 10 + kixa gaps of length k. (b) The path coiresponding to the optimal local alignment generated, superimposed on the hits described in Figure 2. The original BLAST 
program, using the one-hit heuristic with 7- 11, is able to locate three of the five HSPs included in ^s alignment, but only the first and last achieve a score sufficient 
to be reported, (c) The optimal local alignmem, with nominal score 15 and nonnalized score 32.4 bits. In the context of a search of SWISS-PRO'IV(26), release 34 
(21 219 450 residues), using the leghemoglobin sequence (143 residues) as queiy, the £-value is 0.54 if no edge-effect correction (22) is invoked. The original BLAST 
program locates the first and last ungi^iped segments of diis alignment Using sum-statistics with no edge-effect coirection, this combined result has an £-value of 
31 (21^2). On the central lines of the alignment, identities are echoed and substitutions to which tiie BLOSUM-62 manix (18) gives a positive score are indicated 
by a *+• symbol 



reason is that the original BLAST program treats gapped 
alignments implicitly by locating, in many cases, several distinct 
HSPs involving the same database sequence, and calculating a 
statistical assessment of the combined result (21 ,22). This means 
that two or more HSPs with scores well below 38 bits can, in 
combination, rise to statistical significance. If any one of these 
HSPs is missed, so may be the combined result 

The approach taken here allows BLAST to simultaneously 
produce gapped alignments and run significantly faster than 
previously. The central idea is to trigger a gapped extension for 
any HSP tiiat exceeds a moderate score 5g, chosen so that no more 
than about one extension is invoked per 50 database sequences. 
(By equation 2, for a typical-length protein query, 5g should be set 
at -22 bits.) A gapped extension takes much longer to execute 
than an ungapped extension, but by performing very few of them 
the fraction of the total running time they consume can be kept 
relatively low. 

By seeking a single gapped alignment, rather than a collection 
of ungapped ones, only one of the constituent HSPs need be 
located for the combined result to be generated successfully. This 
means fliat we may tolerate a much higher chance of missing any 
single moderately scoring HSR For example, consider a result 
. involving two HSPs, each with the same probability P of being 
missed at the hit-stage of the BLAST algorithm, and suppose that 
we desire to find the combined result with probability at least 



0.95. The original algorithm, needing to find both.HSPs, requires 
2P - < 0.05, or P less than -0.025. In contrast, the new 
algorithm requires only that P'^ < 0.05, and thus can tolerate P as 
high as 0.22. This permits the 7 parameter for the hit-stage of the 
algorithm to be raised substantially while retaining comparable 
sensitivity — ^from 7 = 1 1 to T = 1 3 for the one-hit heuristic. (The 
two-hit heuristic described above lowers Tback to 1 145 As will be 
discussed below, the resulting increase in speed more than 
compensates for the extra time required for the rare gapped 
extension. 

In smnmary, the new gapped BLAST algorithm requires two 
*non-overlapping hits of score at least T, within a distance i4 of one 
another, to^voke an ungapped extension of the second hit If the 
HSP generated has normalized score at least Sg bits, then a gapped 
extension is triggered. The resulting gapped alignment is reported 
only if it has an £- value low enough to be of interest For example, 
in the pairwise comparison of Figure 2, the ungapped extension 
invoked by the hit pair on the left produces an HSP with score 
23.6 bits (calculated using and K^, This is sufficient to trigger 
a gapped extension, which generates an alignment with score 32.4 
bits (calculated using kg and Kg) and £-value of 0.5 (Fig, 3). The 
original BLAST program locates only the first and last imgapped 
segments of this alignment (Fig. 3c), and assigns them ^ 
combined £"value >50 times greater. 
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THE CONSTRUCTION AND STATISTICAL 
EVALUATION OF GAPPED LOCAL ALIGNMENTS 

The standard dynamic prograniming algorithms for pairwise 
sequence alignment perform a fixed amount of computation per 
cell of a path graph, whose dimensions are the lengths of the two 
sequences being compared (23-25). In order to gain speed, 
database search algorithms such as Fasta (2) and an earlier gapped 
version of BLAST (3) sacrifice rigor by confining the dynamic 
programming to a banded section of the full path graph (4), 
chosen to include a region of already identified similarity. One 
problem with this approach is that the optimal gapped alignment 
may stray beyond the confines of the band explored As the width 
of the band is increased to reduce this possibility, the speed 
advantage of the algorithm is vitiated. 

We have accordingly taken a different heuristic approach to 
constmcting gapped local alignments, which is a simple general- 
ization of BLAST'S method for constructing HSPs. The central 
idea is to consider only cells for which the optimal local alignment 
score falls no more than Xg below the best alignment score yet 
found Starting from a single aligned pair of residues, called the 
seed, the dynamic programming proceeds both forward and 
backward through the path graph (Zheng Zhang et a!., manuscript 
in preparation) (Figs 3a and 4). The advantage of this approach 
is that the region of the path graph explored adapts to the 
alignment being constructed. The alignment can wander arbitrari- 
ly many diagonals away from the seed, but the number of cells 
expanded on each row tends to remain limited and may even 
shrink to zero before a boundary of the path graph is encountered 
(Fig. 4). The Xg parameter serves a simUar function to the 
band-width parameter of die earlier heuristic, but the region of the 
path graph it implicitly specifies be explored is in general more 
productively chosen. 

An important elemwit for this heuristic is the intelligent choice 
of a seed Given an HSP whose score is sufficiently high that it 
triggers a gapped extension, how does one choose a residue pair 
to force into alignment? While more sophisticated approaches are 
possible, the simple procedure we have implemented is to locate, 
along the HSP, the length-11 segment with highest alignment 
score, and use its central residue pair as the seed. If the HSP itself 
is shorter than 11, a central residue pair is chosen. For example, 
the first ungapped region in the alignment of Figure 3c constitutes 
the HSP that triggered the alignment The highest-scoring 
length-11 segment of this HSP aligns leghemoglobin residues 
55-65 with P-globin residues 57-67. Thus the alanine residues at 
respective positions 60 and 62 are used as the seed for the gapped 
extension illustrated in Figure 3a. As discussed in the perform- 
ance evaluation section below, this procedure is extremely good 
at selecting seeds that in fact participate in an optimal alignment 

Most gapped extensions are triggered by chance similarities, 
and are therefore likely to be of limited extent, as illustrated in 
Figure 4. The reverse extension in this example explores -2000 
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Figure 4. The path graph region explored by BLAST during a gapped extension 
for the comparison of broad bean leghemoglobin I and the ElB proiein small 
T-anrigen from human adenovirus type 4 (89) (SWISS-PROT accession no. 
PI 0406), The^g dropoff parameter is the nominal score 40, used in conjunction 
with BLOSUM-62 substiwtion scores and 1^ + it gap costs. The 22.7 bit HSP 
that triggers this extension, involving leghemoglobin residues 119-140 and 
adenovirus residues 101-122. is merely a random similarity, and not pan of a 
larger and higher-scoring alignmenL The gapped extension is seeded by the 
alignment of residues 124 and 106. The optimal alignment score through points 
in the path graph drops steadily as one moves beyond tfie triggering HSP, and 
the reverse extension terminates before the beginning of either protein is 
reached. A total of 2766 path gr^h cells are explored, with the reverse 
extension accounting for 2047 of these cells. 



path graph cells, so that a typical two-way gapped extension that 
does not encounter the end of either sequence is expected to 
involve ^000 cells. Because Sg is set so that a gapped extension 
is invoked less than once per 50 database sequences, fewer than 
80 cells need be explored per database sequence. 

The execution time required for a gapped extension is -500 
times that for an ungapped extension. However, by triggering 
gapped extensions in the manner described, while simultaneously 
raising T for the single-hit version of BLAST from 11 to 13, 
approximately one gapped extension is invoked for every 4000 
ungapped extensions avoided. Because the number of ungapped 
extensions is reduced by about two thirds, the total time spent on 
the extension stage of BLAST is cut by well over half. Of course, 
the two-hit strategy described above reduces the time needed for 
the ungapped extensions still further. Once program overhead is 
accounted for, the net speedup is a factor of about three. 

For any alignment actually reported, a gapped extension that 
records *traceback* information (25) needs to be executed. To 
increase BLAST'S accuracy in producing optimal local align- 
ments, these gapped extensions use by default a substantially 
larger parameter than employed during the program *s search 
stage. 



Table 1. Relative times spent by the original and gapped BLAST programs on various algorithmic stages 



Overhead: database 
scanning, output, etc. 



Calculating whether hits Ungapped extensions 

qualify for ung^ped extension 



Gapped extensions 



Original BLAST 



8(8%) 



92 (92%) 
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The times required by various steps of the BLAST algorithm 
vary substantially firom one query and one database to another. 
Table 1 shows typical relative limes spent by the original and the 
gapped BLAST programs on various algorithmic stages. The 
'original BLAST' program is represented, here and below, by a 
variant fomi of blastp version 1.4.9, modified so that it uses the 
same edge-effect correction (22) and background amino acid 
frequencies as the 'gapped BLAST'. The times represent the 
average for three different queries, with the time for the original 
BLAST program normalized in each instance to 100 units. 

More concretely, to search SWISS-PROT (26), release 34 
(59 576 sequences; 21 219 450 residues), with the length-567 
influenza A virus hemagglutinin precursor (27) as query, the 
original BLAST program requires 45.8 s, and the gapped BLAST 
program 15.8 s. This timing experiment, and others referred to 
below, was run on one 200 MHz R 10000 cpu processor of a 
lightly loaded SGI Power Challenge XL computer with 2.5 
Gbytes of RAM. This machine runs the operating system IRK, 
version 6.2, which is an implementation of UNIX. We used the 
standard SGI C compiler, with the -O flag for optimization, to 
compile all versions of the programs. The times reported are the 
user times given by the time command, and are for the better of 
two identical runs. 

A closely related type of gapped extension routine to that used 
here was developed by G. Myers during the evaluation of the 
original BLAST algorithm. It was not included in the publicly 
distributed code primarily because the then current strategy of 
extending every hit decreased the algorithm's speed unduly for die 
relatively small gain in sensitivity realized (1). 

As discussed above, the statistical significance of g^ped 
alignments may be evaluated using the two statistical parameters 
Xg and ^g. The current version of the Fasta program (2) estimates 
these parameters on each run, by analyzing the distribution of 
alignment scores produced by all the sequences in the database. 
BLAST gains speed by producing alignments for only the few 
database sequences likely to be related to the query, and therefore 
does not have the option of estimating and Kg on the fly. 
Instead, it uses estimates of these parameters produced before- 
hand by random simulation (3). A drawback of this approach is 
that the program may not accept an arbitrary scoring system, for 
which no simulation has been performed, and still produce 
accurate estimates of statistical significance. The original BLAST 
programs, in contrast, because they dealt only with ungapped 
local alignments, could derive and from theory for any 
scoring matrix (8,9). 

ITERATED APPLICATION OF BLAST TO 
POSITION-SPECIFIC SCORE MATRICES 

Database searches using position-specific score matrices, also 
called profiles or motifs, often are much better able to detect weak 
relationships than are database searches that use a simple 
sequence as query (28-38). Employing these methods, however, 
frequently has involved the use of several different programs and 
a fair degree of expertise. Accordingly, to render the power of 
motif searches more readily available, we have written a 
procedure to construct a position-specific score matrix automati- 
cally from the output of a BLAST run, and modified BLAST to 
operate using such a matrix in the place of a simple query. The 
resulting PSI-BLAST program often is substantially more 
sensitive than the corresponding BLAST program, but for each 



iteration takes litde more than the same time to run. In related 
work, Henikoff and Henikoff (39) have described how, short of 
modifying BLAST so that it may operate on a position-specific 
score matrix, a single artificial sequence that approximates such 
a matrix may be used as a query with the original BLAST 
programs. 

The construction of a position-specific score matrix is a 
multi-stage process, and at each stage a choice must be made 
among a number of alternative routes. We have been guided by 
the goals of automatic operation, speed of execution, and general 
simplicity. The issues discussed below are: (i) general architec- 
ture of the score matrix; (ii) construction of the multiple 
alignment from which the matrix is derived; (iii) weights for 
sequences within the multiple alignment, and evaluation of the 
effective number of independent observations it constitutes; 
(iv) estimation of target frequencies, and the construction of 
matrix scores; (v) applying bL\ST to a position-specific matrix, 
and the statistical evaluation of search results. We do not claim 
our current implementation is optimal, and it is likely that over 
time some of its details will change. 

Score matrix architecture 

The alignment of a simple sequence with a pattern embodied by 
a position-sp)ecific score matrix is almost completely analogous 
to the alignment of two simple sequences. The only real 
difference is that the score for aligning a letter with a pattern 
position is given by the matrix itself, rather than with reference to 
a substitution matrix. For proteins, a query of length L and a 
substimtion matrix of dimension 20 x 20 are^ replaced by a 
position-specific matrix of dimension L x 20. Position-specific 
gap costs may be defined as well (34,40). As with pairwise 
sequence comparison, one may choose among finding the best 
global alignment of the matrix and the simple sequence (23), 
finding the best alignment of the complete matrix with a segment 
of the sequence (41), and finding the best local alignment of the 
matrix and sequence (24). 

Position-specific protein score matrices draw their power from 
two sources. The first is improved estimation of the probabilities 
with which amino acids occur at various pattern positions, leading 
to a more sensitive scoring system. The second is relatively 
precise definition of the boimdaries of important motifs. By 
demanding the complete aiigrunent of one or more iHotifs, rather 
than seeking an arbitrary local alignment, the size of the searcii 
space may be greatly reduced, thereby lowering the level of 
random noise. Unfortunately, there are many obstacles to 
automating well the delineation of a set of motifs from the output 
of a database search. The query sequence may contain a variety 
of differejjt domains, and share different subsets of them witfi 
different proteins in the database. Furthermore, defining the 
proper extent of even a single motif may be challenging (42). 

Accordingly, we have chosen to forgo the potential advantages 
of restricting the length of our derived matrices, and ihcn 
demanding that they be completely aligned with segmenvs of 
database sequences (41). Instead, each matrix we construct has 
length precisely equal to that of the original query sequence. 
When searching the database with such a matrix, we seek local 
alignments, in full analogy to those sought by BLAST when used 
for straightforward sequence-sequence comparison. Finally, we 
do not attempt to derive position-specific gap scores for use with 
our position-specific substitution scores. Instead, in each iteration 
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3 Accession 



Alignment 



i. 



1 




Bistidl&e triad protalB 15 vflktslsfalvnrkpwpchvlvcplrpvsrphdlrpdev^dlp 59 

♦ ++TB ALV t P l** P V-fR +L -M. DL 
OridylyltreaBferase 213 IVVBTEHMIALVPYVXlHPFBTLLLPKIHVKIttTELSDEOSimui 257 



R&Btidine triad protein 
Dridylyl transferase 



60 OTTORVCTWBICHFHCT-SI.TPSMODCPEACQTVKH--VHVHVLP 101 
256 VZLKRLTTKyDNLFETSPPYSHGPHAAPFNSEDMSHWQLHJiBFTP 302 



Biattdine triad protein 102 R--iaiCDPBIlKDSirZELOKHDKEDrPASWRSEEEMAAEAAALRV 144 

♦ YE L -M- + "M- XE AA R+ 

Orldylyltranaferase 303 PLLRSATVRKTMVCyEMLGEN ORDLTAEQAAERL 336 



BiBtidine triad prorein 
Phosphorylase 

Histidine triad protein 
Phosphorylase 



25 LVNRKPWPGHVLVCPLRPVERFHDLRPDCVADLPQTTORVGT 67 

L+N+ PV+PCH L"^ + L p +^ T 
91 LLNKFPVIPGHTLLVTNEYQHQTDALTPTDL LTAYKLLC 129 

66 WEKBPHCTSLTrSMODGPEACQTVKHVHVHVL--PRKAGDF 107 
-M- CP +0 ♦'f H H+ *L P K F 

130 ALDNEESDKRHMVFyKSGPASGSSLDHKHt.0JLOMPEKFVTF 171 



Sp^Vto S^S^^^^^^ generated by PSI-BLAST when the human ftagilc histidine triad (HIT) protein (61) (SWISS-PROT accession no. P49789) is 

STX^rlc^^"^^^^^ P^''" alignments have £-value SO.OI . and are identified in SWISS-PROT as belonging to the HIT family. TOck bars within 

nlX^n^^I^^^T ^ ^^^^dedpomons of the multiple angnmemaie used (b)AlocalaHgnmem of the hum and 

^^nfluen^e galactose- 1 -phosphate undylylinmsferase (63) (SWISS-PROT accession no. P31764). In its first position-spe^ific iteration, PSl-BLAST riv« Ss 

no. P16550). In its second position-specific iteration, PSI-BLAST gives this aUgnment a score of 43,4 bitsTcoiresponding to an £-value of 2 x 1(H 



of PSI-BLAST, we employ the same gap scores that are used in 
the first, simple BLAST run. Our reasons are that there is no good 
theory for deriving gap costs from a multiple alignment and that, 
as will be discussed below, by eschewing position-specific gap 
costs we can make a reasonable estimate of the statistical 
significance of the resulting -local alignments. 

Multiple alignment construction 

To produce a multiple alignment from the BLAST output, we 
simply collect all database sequence segments that have been 
aligned to the query with £-value below a threshold, by default 
set to 0.01. The query is used as a master, or template, for 
constructing a multiple alignment M. Any row (i.e., database 
sequence segment) identical to the query segment with which it 
aligns is purged, and only one copy is retained of any rows that 
are >98% identical to one another. Pairwise alignment columns 
that involve gap characters inserted into the query are simply 
Ignored, so that M has exactly the same length as the query. 
Because we are dealing with local alignments, the columns of M 
may involve varying numbers of sequences, and many columns 
may include nothing but the query. We make no attempt to 
improve M by comparing database sequences with one another, 
or by any other true multiple alignment procedure. • 

As will be discussed, the matrix scores constructed for a given 
alignment column should depend not only upon the residues 
appearing there, but upon those in other columns as well. To make 
this dependency easy to formulate, however, we need to prune our 
i^w multiple alignment M to a simpler ^reduced* one. This 
Pnaning is done independentiy for each column, so the reduced 
multiple alignment Mc wiU in general vary from one column C 
^0 the next To construct Mc, we first specify the set R of 
sequences it includes to be exactly those that contribute a residue 
to column C. We then define the columns of Mc to be just those 
columns of M in which all the sequences ofR are represented. By 
construction, the reduced multiple alignment Mc has residues or 



gap characters in every row and colunm (Fig. 5a), and is therefore 
amenable to the various manipulations described below. 

Sequence weights 

When constructing a score matrix from a multiple alignment, it 
is a mistake to give all sequences of the alignment equal weight. 
A large set of closely related sequences carries little more 
information than a single member, but its size alone may allow it 
easily to 'outvote* a small number of more divergent sequences. 
One way past this difficulty is to assign weights to the various 
sequences, with those having many close relatives receiving 
smaller weight The many sequence weighting methods that have 
been proposed (43-5 1) often produce roughly equivalent results. 
Because of its speed and simplicity, we have implemented a 
modified version of the sequence weighting method of Henikofif 
and Henikoff (47). Gap characters are treated as a 21st distinct 
character, and any columns consisting of identical residues are 
ignored in calculating weights. In speaking of a column's 
obsCTved residue frequencies /j, we shall henceforth mean its 
weighted rather its raw frequencies. 

In constructing matrix scores, not only a column's observed 
residue frequencies are important, but also the effective number 
of independent observations it constitutes: a column consisting of 
a single valine and a single isoleucine carries different informa- 
tion than one consisting of five independentiy occurring instances 
of each. Accordingly, we need to estimate tiie relative niunber A^c 
of independent observations constituted by the aligrunentMc. A 
simple count of the number of sequences in Mq is a poor measure, 
for 10 identical sequences imply fewer independent observations 
than do 10 divergent ones. We thus propose as a simple first 
estimate for Nq die mean number of different residue types, 
including gap characters, observed in the various columns of Mc- 
This estimate is clearly not ideal, as it saturates at 21 no matter 
how many independent sequences are contained in Mc- However, 
for die data we are likely to encounter. Nq is typically much 
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smaller than 21 , and therefore perhaps a good enough approxima- 
. tion for our puiposes. As will be seen, it is not the absolute value 
of A^c is important, but rather its relative value from one 
column to another. Nq is essentially the same measure of 
alignment variability as that proposed by Henikoff and Henikoff 
(52) for use in a different manner. 

Target frequency estimation 

Given a multiple alignment, many methods for generating score 
matrices have been advanced (28-37,42,52-54). The prescrip- 
tion with perhaps the best theoretical foundation is that the scores 
for a specific pattern position be of the form log {Qi/PO, where Qi 
is the estimated probability for residue i to be found in that column 
(29,30,32,33,36,37,42,52-54). This leaves open the question 
how best to estimate the Q\. 

Given a multiple alignment involving a large number of 
independent sequences, the estimate of Qi for a specific column 
should converge simply to the observed frequency of residue i in 
diat column. However, in addition to the sequence weighting 
issues discussed above, factors that complicate estimating the Q\ 
include small sample size (30), and prior knowledge of relation- 
ships among the residues (1 6,37,53). Various studies suggest that 
the best currently available method for estimating the Qi is that of 
Dirichlet mixtures (52-56). However, because it often perfomis 
nearly as well (52), and due to its relative simplicity, we have 
implemented the data-dependent pseudocount method intro- 
duced by Tatusov et al. (31). liiis method uses Ae prior 
knowledge of amino acid relationships embodied in the substitu- 
tion matrix Si^ to generate residue pseudocount frequencies gi, 
which are averaged with the observed frequencies^ to estimate 
the Qi, 

Specifically, for a given column C, we construct pseudocount 
frequencies gi using the formula: 

J ' 

where the are the target frequencies implicit in the substitution 
matrix, and given by equation 3. Intuitively, those residues 
favored by the substimtion matrix to align with the residues 
actually observed receive high pseudocount frequencies. We then 
estimate Qi by: 

where a and p are the relative weights given to observed and 

pseudocount residue frequencies. So that the scores we constmct 
will reduce to in colunms where nothing has been aligned to the 
query sequence, we let a = Nc - 1- [3 remains an arbitrary 
pseudocount parameter, the larger its value, the greater the 
emphasis given to prior knowledge of residue relationships vis a 
vis observed residue frequencies. We have found empirically that, 
in conjunction with our method for calculating a, a reasonably 
good setting for P is 10. 

BLAST applied to position-specific score matrices 

The initial step of the BLAST algorithm is the construction of a 
list of words that align to query words with score at least T. Only 
minor modifications to the code are necessary for this step to be 



performed on a query consisting of a position-specific matrix 
rather than a simple sequence. The same holds for the ungapped 
and gapped extension, steps of BLAST One important issue is 
whether key parameters such as TandXg, used at various heuristic 
stages of the algorithm and tuned to simple sequence comparison, 
can be applied unchanged to position-specific matrices without 
degrading unduly either the speed or sensitivity of database 
searches. We approach this problem by ensuring that the scale 
of the matrix scores produced internally by PSI-BLAST corre- 
sponds to that of the substimtion matrix Siy In other words, we 
calculate the scores for a column of the matrix as []n(Qi/Fi)]/7^, 
There is no analytic . theory with which to estimate the statistical 
significance of a gapped alignment of a position-specific score 
matrix and a simple sequence. However, one may hypothesize 
that for a score matrix constructed to the same scale as 5"ij, a given 
set of gap costs should produce the same gapped alignment scale 
parameter as for s\y This would be convenient, because then 
PSI-BLAST could estimate statistical significance without 
expending after each iteration the substantial time required to 
estimate ^ and Kg by random simulation. To test this hypothesis, 
we performed a number of statistical tests on PSI-BLAST 
generated score matrices, scaled to have Xa - 0.3176, the value 
applicable to previously published BLOSlJM-62 simulations (3). 

First, we searched S WISS-PROT using as query the lengtii-567 
influenza A virus hemagglutinin precursor (27), and captured the 
score matrix constructed by PSI-BLAST from the 128 local 
alignments witii £-value < 0.01. We then compared this matrix to 
10 000 random sequences of length 567, generated using the 
background amino acid frequencies of Robinson and Robinson 
(20). A gap of length k was charged a cost of 10 + k. Counts of 
the optimal local alignment scores, calculated using an appropri- 
ately modified version of the Smitii-Waterman algorithm (24), 
are plotted in Figure 6. Also shown is the best fitting extreme 
value distribution (3,15) which, using the edge-effect correction 
described by Altschul and Gish (3), has statistical parameters 7^ 
= 0.25 1 and Kg = 0.03 1 . It is apparent that the distribution fits :he 
random trial reasonably well; a goodness-of-fit test with 34 
degrees of freedom has value 41.8, which is lower than one would 
expect 20% of the time even were the theory precisely valid. This 
supports the idea that die statistical theory described above 
applies to local alignments of position-specific score matrices and 
simple sequences. Furthermore, the estimate Xg = 0.251 ± 0.003 
agrees to within experimental error with the value 0.255 
previously published for these gap costs (3). Similar agreement 
was obtained with a number of other protein sequen^s as initial 
query (results not shown), and in all cases the much less important 
Kg parameter could be estimated accurately as well. In general, 
vines of Xg for the comparison of position-specific score 
matrices with simple sequences appear to differ by <2% from the 
values for simple pairwise sequence comparison. Using these 
«precomputed values for Xg should thus entail an error of less than 
one bit for PSI-BLAST scores <:5G bits, corresponding to a factor 
of less than two in the estimation of statistical significance. 

PERFORMANCE EVALUATION 

To test more directiy the statistics used by PSI-BLAST, we 
compared query sequences from 1 1 large and well characterized 
protein families to tiie SWISS-PROT database, and then ran the 
position-specific score matrices generated against a shuffled 
version of the same database. For each query, we recorded the 
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1000 



Counts 




Optimal local alignment 



Figure 6. The disnibution of optunal local alignment scores from the 
comparison of a position-specific score matrix with 10 000 random protein 
sequences. The score matrix was constructed by PSI-BLAST from the 128 local 
alignments with £-value ^.01 found in a search of SWISS-PROT using as 
query the length-567 influenza A virus hemagglutinin precursor (27) (SWISS- 
PROT accession no. P03435). The random sequences^ each of length 567, were 
generated using the amino acid frequencies of Robinson and Robinson (20). 
Optimal local alignment scores were calculated using the position-specific 
matrix in conjunction with 10 + ^.gap costs. The extreme value distribuDon that 
best fits the data (3,15) is plotted. A goodness-of-fit test with 34 degrees of 
freedom has value 41.8, corresponding to a /'-value of 0.20. 

lowest £-vaIue found, as well as the number of shuffled sequences 
yielding £-values <1 and 10. For comparison, we performed the 



identical shuffled-database test on the gapped and original 
versions of BLAST. To reduce the probability that high-scoring 
alignments were missed due to the heuristic nature of the 
algorithms, we perfonned these tests with T = 9 rather than the 
default value of 11. The results are given in Table 2. For the 11 
queries, the median of the low PSI-BLAST £-values was 0.87, 
which corresponds to a median P-value of 0.58 (8,9). The mean 
numbers of shuftled database sequences with £-values <1 and 10 
were 1 .0 and 8.7, respectively, within 20% of the expected values 
of 1 .0 and 10.0. The equivalent tests for the ungapped and gapped 
versions of BLAST also yielded results thai diverged from theory 
by<50%. 

The ability to estimate with reasonable accuracy the signifi- 
cance of gapped local matrix-sequence alignments permits us to 
automate the construction of position-specific score matrices 
during multiple iterations of the PSI-BLAST program. After each 
iteration, we generate a new multiple alignment simply by 
collecting those alignments with £-value lower than a defined 
threshold. An interactive version of PSI-BLAST allows the user 
to override either the inclusion or exclusion of specific local 
alignments. Once a given database sequence has been used in the 
generation of a position-specific score matrix, low £-values for 
this sequence are virmally guaranteed in future iterations, for the 
sequence is to a certain extent being compared with itself. The 
biological relevance of PSI-BLAST output thus depends criti- 
cally on avoiding the inappropriate inclusion of sequences in the 
multiple alignment constructed. Specifically, the utility of the 
score matrix produced is immediately vitiated by the inclusion of 
any alignment involving a region of highly biased amino acid 
composition (57,58). 



Table 2. The comparison of various query sequences with a shuffled version of SWISS-PROT 



Protein family 


SWISS-PROT 


Original BLAST 




Gapped BLAST 




PSI-BLAST 






accession no. 


Low 


No. of 


seqs 


Low 


No. of seqs 


Low 


No. of seqs 




of query 


£-value 


with £-va]ue 


£-value 


with £-value 


£-value 


with £-valuc 








^1 


^10 




SI 


<10 




SI 


<10 


Serine protease 


P00762 


0.86 


1 


7 


3.0 


0 


4 


0.94 


1 


8 


Serine protease inhibitor 


P01008 


3.9 


0 


4 


0.078 


1 


9 


1.5 


0 


9 


Ras 


poini 


3.4 


0 


8 


3.4 


0 


7 


1.1 


0 


9 


Globin 


P02232 


2.4 


0 


7 


2.8 


0 


5 


8.2 


0 


2 


Hemagglutinin 


P03435 


0.11 


2 


11 


0.46 


3 


16 


0.87 


1 


8 


Interferon a 


P05013 


2.4 


0 


6 


0.27 


2 


4 


0.11 


2 


11 


Alcohol dehydrogenase 


P07327 


1.5 


0 


2 




1 


5 


1.5 


0 


9 


Histocompatibility antigen 


P10318 


0.91 


1 


7 


0.13 


1 


7 


0.0031 


2 


6 


Cytochrome P450 


P10635 


0.84 


2 


5 


8.5 


0 


3 


0.46 


1 


15 


Glutathione transferase 


P14942 


1.0 


1 


10 


3.3 


0 


3 


0.30 


2 


9 


H"*--transporting ATP synthase 


P20705 


0.012 


1 


8 


0.26 


2 


14 


0.79 


2 


10 


Average (median or mean) 




LO 


0.7 


6.8 


0.80 


0.9 


7.0 


0.87 


1.0 


8.7 



original and gapped BLAST comparisons use BLOSUM-62 substitution scores (18). All three programs use threshold T parameter set to 9, but the gapped 
BLAST and PSI-BLAST programs use the two-hit method to trigger ungapped extensions. The original BLAST program has the X dropoff parameter set to nominal 
score 23. Hie gapped BLAST and PSI-BLAST comparisons charge gaps of length k a cost of 1 0 + Jk. They have Xu set to 1 6. and X^ set to 40 for the database search 
stage and to 67 for the output stage of the algorithms. Gapped alignments are triggered by a score corresponding to -22 bits. For PSI-BLAST, the query is first com- 
pared to the SWISS-PROT database, and the position-specific score matrix generated is then compared to a shuffled version of SWISS-PROT. The median is used 
for the average of the low £-values, and the mean otherwise. 
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Table 3. The number of SWISS-PROT sequences yielding aiignmenu with £-value ^0.0 1 , and reiaii ve running limes, for Smith-Waterman and various versions 
of BLAST 



Protein familv 


Query 


dmi ui^ wacerman 




uapped cLAb I 


DOT DT A PT* 

rSi-BLAST 


Serine protease 


rUU702 


lie 

275 


273 


275 


286 


Serine protease inhibitor 


P01008 


108 


105 


108 


111 


Ras 


ponii 


255 


249 


252 


375 


Globin 


P02232 


28 


26 


28 


623 


Hemagglutinin 


P03435 


128 


114 


128 


130 


Interferon a 


P05013 


53 


53 


53 


53 


Alcohol dehydrogenase 


P07327 


138 


128 


137 


160 


Histocompatibility antigen 


P10318 


262 


241 


261 


338 


Cytochrome P450 


PI 0635 


211 


197 


211 


224 


Glutathione transferase 


PI 4942 


83 


79 


81 


142 


H^-transporting ATP synthase 


P20705 


198 


191 


197 


207 


Normalized running time 




36 


1.0 


0.34 


0.87 



To score and evaluate the significance of the alignments found, the original BLAST program uses BLOSUM-62 substitution scores (18) and sum-statistics (21,22), 
The Smith-Waterman and gapped BLAST programs use BLOSUM-62 subsUtuiion scores, 1 0 + Jk gap costs/and the statistics of equations 1 and 2, in conjunction 
with the experimentally determined parameters = 0.255 and ATg = 0.035 (3). PSI-BLAST uses the same gap costs and Xg, but applied to the position-specific score 
matrix constructed from the output of the gapped BLAST mn. Only one PSI-BLAST iteration is executed. All three BLAST programs use the same parameter settings 
as in Table 2, except that Tis set to 1 1. Normalized running times are the mean ratio of program running time to that for the original BLAST. The time for PSI-BLAST 
includes the time for the initial BLAST search. 



To compare the performance of the new gapped version of 
BLAST and its PSI-BLAST extension to that of the Smith- 
Waterman algorithm (24) and the original ungapped BLAST 
algorithm, we employed the same 11 query sequences that were 
used above to investigate the accuracy of PSI-BLAST statistics. 
Because, as shown, these statistics are quite accurate, we may use 
the number of statistically significant sequences found in a 
database search as a reasonable measure of algorithm sensitivity. 
We employed the ssearch program, version 2.0u54, from the 
Fasta package (2) as our implementation of the Smith-Waiennan 
algorithm. Using each of the 11 queries, we searched SWISS- 
PROT with each of die four programs. We show in Table 3 the 
numbers of sequences found with E value <0.01, as well as the 
average ratio of running time to that for the original BLAST 
program. Based upon SWISS-PROT annotation, all sequences 
recorded in Table 3 appear to be true family members, with the 
exception of one of iht lowest-scoring alignments found by 
Smith-Waterman when applied to the histocompatibility antigen 
query, and the lowest-scoring alignment found by the original 
BLAST applied to the hemagglutinin query. While some 
alignments involve hypothetical proteins, the pattern of con- 
served residues in all such cases suggests a true positive. 

As can be seen, the gapped BLAST program runs on average 
three times faster then the original, and in all but one case 
examined finds a greater number of statistically significant 
alignments. It runs >100 times faster than Smith-Waterman, but 
for the combined 11 queries misses only eight of the 1739 
significant similarities found by the rigorous algorithm. Of these 
eight, only one has an £-value <0.001, and another ^qjpears to be 
a random as opposed to a biologically meaningful similarity. The 
scores produced by gapped BLAST for the 1731 similarities it 
fmds^ differ from those produced by the Smith-Waterman 
algorithm in only two instances. The discrepancy arises in both 
cases from an Xg parameter that is too low rather than from an 
incorrect choice of seed. Thus despite its simplicity, the 
seed-selection heuristic is extremely accurate. 



A search that includes a single PSI-BLAST iteration still runs 
faster than the original BLAST, and 40 times faster than 
Smith-Waterman, but can in many cases be much more sensitive. 
It finds every true positive returned by Smith-Waterman, but 
frequendy many others as well. Here only a single PSI-BLAST 
iteration has been considered but, as will be seen below, multiple 
iterations can yield even better results. Furthermore, we have 
found PSI-BLAST to perform better on searches of the non- 
reduridant protein sequence database maintained by the NCBI 
(59) than on searches of SWISS-PROT. because of the greater 
number of significant similarities that arc found by the initial 
BLAST run. 

For tiie particular examples in Table 3, the PSI-BLAST 
iteration takes noticeably longer than the gapped BLAST 
iteration, due primarily to the time needed to construct the 
position-specific score matrix from the large number of signifi- 
cant local alignments found by BLAST. For queries that return a 
small number of significant alignments, each F^I-BLAST 
iteration requires more nearly the same time as BLAST 

PSI-BLAST EXAMPLES 

In many instances, PSI-BLAST is able automatically to uncover 
biologically* interesting similarities that elude simple database 
searches. Multiple iterations of PSI-BLAST are sometimes required 
to recognize the more disiantiy related protein family members. We 
here consider two representative cases in greater detail 

HIT proteins 

Holm and Sander (60) describe how a comparison of three-dimen- 
sional structures is able to identify significant similarity between 
histidine triad (HIT) proteins and galactose- 1 -phosphate uridylyl- 
transferase (GalT) proteins. Indeed, using the human HIT protein 
(61) as query, a BLAST search of SWISS-PROT reveals hits.with 
£-value <0.01 only to other HET proteins (Fig. 5a). An alignment to 
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^ the rat GalT protein (62) has tfie only marginally significant E-value 
of 0.012. A PSI-BLAST search, using the score matrix generated 

I from the six alignments illustrated in Figure 5a, can immediatdy 

1 cement confidence in the biological idevance of this similarity. The 

I £^ueofthesinnilarity withratGalTdropsto2x 10"^, andan 

I alignment to Haemophilus influenzae GalT (63) (Kg. 5b) receives 

I the even more significant £-value of 4 x 10"^. These similarities, of 

r- course, are uncovered using no structural infommon. In addition, 

I on the next iteration, PSI-BLAST finds a strongly significant 

I alignment (Hg. 5c; £Walue 2 x 10^ to yeast 5\5'":P1 J^4-tetia- 

I phosphate phosphoiylase I (64), for which no structure is available. 

BRCT proteins 

I Proteins containing one or multiple copies of the BRCT domain 
I form a superfamily many of whose members are involved in DNA 
I damage-responsive cell cycle checkpoints (65-67). While detailed 
I analysis is needed to delineate completely Ais diverse set of 
I proteins, PSI-BLAST is able to automatically identify most of the 
superfamily. We used the C-terminal 215 residues of human 
BRCAl (68), which includes two BRCT domains (65), as the 
initial queiy for a search of NCBFs non-redundant protein 
sequence database. Using the default cutoff £-value of 0.01, the 
initial BLAST search recognized as significant only alignments to 
other BRCAl sequences, and the previously described BRCT 
protein BARD (69) (Table 4). Subsequent PSI-BLAST iterations, 
however, retrieved all the proteins recorded in Table 4; additional 
close homologs are omitted from the table. Ahnost all the BRCT 
proteins described by Bork et al. (66) were recognized. Not found 
were the retinoblastoma family, whose putative BRCT domain is 
particularly diveigent, worm R13A5.13, which was not in the 
database searched, and human DNA-ligase m. PSI-BLAST did 
report yeast RAD9 and YGR103w, the Kluyveromyces lactis 
RAPl homolog, womi ZK675.2, and various temiinal deoxynu- 
cleotidyltransferases and poly(ADP-ribose) polymerases, but all 
with £-values >0.01 (Table 4). Detailed examination of the 
alignments produced suggests that the only likely false positives 
involved a tiypanosome EST (70) and the Methanococcus 
jannaschii mutT protein (71), the latter despite its involvement in 
DNA repair (Table 4). 

Seven recent additions to the protein databases are reported here 
as members of the BRCT superfanuly (Table 4). (i) Arabidopsis 
T10M13.12 (72), is the first plant protein observed to contain 
BRCT domains, (ii) KIAA0259 (73) is a large human protein of 
unknown function with eight BRCT domains, tiie greatest number 
so far observed within a single protein, (iii) Tl 3F2.3 (74) is a worm 
protein with a 500-residue low-complexity (57) N-terminus. 
(IV) SPAC6G9.12 (75) is a fission yeast protein strongly similar to 
the previously recognized (66) yeast BRCT protein L8543. 1 8 (76). 

(v) C36A4.8 (74) is a worm protein whose C-teiminus contains a 
single BRCT domain, and whose N-terminus, containing a RING 
finger domain, is strongly similar to tiiat of BRCAl. The similar 
organization to BRCAl makes this protein of particular interest 

(vi) Synechocyms sp. D90904 (77) is tiie first bacterial BRCT 
protein tiiat is not a bacterial ligase. While it failed to pass tiie cutoff 
^'-value of 0.01, its C-terminal BRCT domain is very similar to that 
of several bacterial ligases, which presumably led to its incon^t 
^sification as such in the databases. Most of die protein 
^terminal to its BRCT domain consists of a .coiled-coil domain. 
The actual Synechocystis sp, DNA ligase (77.78) is found by 
^SI-BLAST on the 13th iteration, witii an £-value of 0.002. 



H.sapUm BRCAl Q 

A.thaUaia T10M13.12 M . aBCTi 



H, sapiens KiAA0259 
Cekgans T13FZ3 
S,pombe SPAC6G9.12 




Celegans C36A4.B Q £3 
Synechocystis sp. D90904 ^ 

H. sapiens Pescadillo gg) JOO 

Figure 7. The location of BRCT domains within human BRCAl (68), 
Arabidopsis T10M13.12 (72), human KIAA0259 a3), worm T13F23 (74). 
fission yeast SRAC6G9.12 (75). worai C36A4.8 (74), Synechocystis sp. 
D9()904 (77) and human Pescadillo (79). BRCL\1 and C36A4.8 each have, in 
addition, an N-temiinal RING finger domaia The near ideniiTy to other worm 
sequences of a shon region directly preceding the BRCT domain of C36A4.8 
suggests tiie possibility that this protein has been nusasscmbied. 

(vii) Pescadillo is a human protein whose zebrafish ortholog is 
essential for embryonic development (79), and whose yeast 
ortholog YGR103W (80) has been previously recognized as a 
BRCT protein (66,67). It failed to pass cutoff £-value of 0.01, 
but appeared with near-significant £-values in PSI-BLAST output 
from tiie 5th iteration onward. The approximate positions of the 
BRCT domains within BRCAl and the seven newly identified 
BRCT proteins are illustrated in Figure 7. 

DISCUSSION AND CONCLUSION 

In addition to the major algorithmic changes described above, we 
have modified an aspect of the original BLAST program's o\xcpal 
routine that on occasion caused important similarities to be 
overlooked. When a very large number of statistically significant 
alignments was found, BLAST would typically report only the top 
scoring 500. These alignments, however, might all involve one 
domain of the queiy that occurred frequently within the database. 
Interesting but weaker relationships to other regions of the query 
might simply be forced off the bottom of the^lisL Accordingly, 
following the general idea of Sonnhammer and Durbin (81), we 
have limited the number of alignments reported that involve each 
region of the queiy, but set no ov«^ upper limit. 

The BLAST programs are unlikely to remain static, and there 
are many possible avenues for future improvement We discuss 
three of them briefly here. 

Gap costs 

Gapped alignments may be constmcted using a variety of different 
types of gap cost Because a single mutational event may insert or 
delete a large number of residues, it has been argued that long gaps 
should not cost much more than short ones, and cffim gap costs, 
which assess a score -(a + bk) for a gap of length k (82-85), have 
become the most widely used. A generalization of these costs has 
been proposed, that aUows a gap to involve residues in both 
sequQices rather than just one (86). Specifically, a gap in which k . 
residues are inserted or deleted and j pairs of residues are left 
unaligned receives the score -<<2 + + cf). The algorithm necessary 
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for using such costs is only a minor variant on that for traditional 
affine gsp costs. In many cases, the new gap costs generate local 
alignments that are bodi more accurate and more statisdcally 
significant (86). These costs are potentially of particular value for 
use with PSI-BLAST, because by in?)Osing alignment only where 



it is justified, they may lead to the construction of more sensitive 
position-specific score matrices. Whedier it is desirable to use 
generalized aSBne g^ costs as the default for geneial purpose 
daiabase searches awaits detailed empirical study. . 



Table 4. PSI-BLAST protein database search results using the C-ierrainus of BRCAl as query 



Protein 



Species 



GenBank ID number PSI-BLAST iteration £-value 



BARD 


Homo sapiens 


1710175 


0 


2e-06 


T10M13.12* 


Arabidopsis thaliana 


2104545 


1 


4C-06 


F26D2.bb 


Caenorhabditis elegans 


1914176 


1 


4e-04 


KIAA0259« 


H.sapiens 


1665785 


1 


0.001 


F37D6.1 


C. elegans 


1418521 


2 


4e-06 


C19G10.07 


Schizosaccharomyces pombe 


1723501 


2 


6c-05 


KIAA0170 


H^apiens 


1136400 


2 


0.002 


53BP1 


H.sapiens 


488592 


2 


0.008 


T13F2.3* 


Celegans 


1667334 


3 


2C-07 


K04C2.4 


Celegems 


470351 


3 


3e-07 


T19E10.1 


. Celegans 


1067065 


4 


7e-04 


Rad4/Cut5 


S.pombe 


730470 


4 


0.002 


REVl 


Saccharomyces cerevisiae 


. 132409 


4 


0.003 


ECT2 


Mus musculus 


423597 


5 


le-04 


XRCCl 


Mmuscuius 


,627867 


5 


6e-04 


Crb2 


o.poTnoe 


1 Add 7*7 


c 


\J.\J\JXm 


P API 


u tCcrcvisuie 


1 / J J JO 


K 
tj 


0.006 


ICJfco 1 u^lr 


Tryponosoftta cruzi 




D 


u.w/i. 


UrDLl 


U . CcicVl 4 lOc 


1352999 


g 


0.001 


T QCA'X 1 Q 


S.cerevisioe 


lU/OU/ J 




0.010 




S^pOTTtbe 




7 


4e-04 




S .cerevisiae 




7 


0.005 


VTJD 1 Ci4»r 

inKiD4W 


Sxerevisiae 




7 




C36A4.8* 


Celegans 


1657667 


7 


0.010 


UNE452 


S.cerevisiae 


1151000 


8 


8e-04 


DNA Ugase IV 


H^apiens 


1706482 


8 


0.008 


CDC9 


Candida albicans 


1706483 


9. 


0.006 


DNA ligase 


Thermus scotoductus 


1352293 


10 


0.010 


GNFl 


Drosophila melanogaster 


544404 


11 


0.004 


inutT^ 


Mjannaschii 


2129134 


15 


0.008 


RAD9 


Sxerevisiae 


131817 


7 


.0.74 


RAPl homolog 


KMctis 


422087 


9 


' 0.21 


ZK675.2 


Celegans 


599712 


13 


3.5 


D90904* 


Synechocystis sp. 


1652299 


15 


0.17 


TDT . 


Mus domestica 


2149634 


15 


o:46 


yGR103w 


S.cerevisiae 


172»693. 


16 


0.017 


Pescadillo* 


H^apiens 


2194203 


16 . 


0.017 


PPOL 


Sarcophaga peregrina 


1709741 


16 


0.060 



Iteration zero refers to the initial BLAST nm, using the 215 C-terminal residues of BRCAl (68) (SWISS-PROT accession no. 
P38398) as query. Subsequent PSI-BLAST iterations use derived position-specific score matrices in place of the query. The score 
matrix for iteration / + 1 is constructed from aligmnents achieving an £-value ^0.01 for iteration /, For each protein, the £-value is 
that returned during the PSI-BLAST iteration indicated, and precedes the protein's use for score matrix construction. Only one repre- 
sentative is listed for families of closely related proteins. On its 16th iteration PSI-BLAST uncovered no new proteins with £-value 
<0.01 , and therefore ceased iteration. At the end of the table are sho>^'n BRCT proteins remmed by PSI-BLAST with £-value >0.01 
but <10, listed for the iteration in which they achieved their lowest £- value. 
^Recent additions to the database, first identified as BRCT proteins here. 

^e Celegans F26D2.b protein (74) while a recent addition to the databases, is a close homolog of the previously recognized (66,67) 
family of C.c/cganj BRCT proteins containing, for example, F37A4.4 (90), 

^e trypanosome EST (70) and the Mjannaschii mutT protein (7 1) are the only likely false positives. 
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Position-specific score matrices as input to PSI-BLAST 

FSI-BLAST performs three distinct operations: it constmcts a 
multiple alignment from BLAST output data; it processes this 
alignment into a position-specific score matrix; and it uses this 
matrix to search the database. A researcher may wish, however, to 
bypass iStit first two of these op^ations, and provide a score matrix 
as query directly to PSI-BLAST. Hie central difficulty is retaining 
the abUity to calculate reliable statistics; as described above, 
PSI-BLAST imposes strict scaling rules on the matrices it generates, 
pennitdng the use of precomputed Xg to assess significance. Tbree 
possible routes are open, (i) (Dne may permit the spedficaiion of 20 
target frequencies 2i for each position rather than 20 scores. These 
can then be converted internally to Ipg-odds scores with the 
appropriate scale for precomputed . statistical parameters to ^ly. 

(ii) One may estimate by random simulation the statistical para- 
meters for an input score matrix (3). An advantage is the 
e^licability to a greater range of scoring systems, including the 
possible use of position-specific gap costs. A disadvantage is that 
obtaining reasonably accurate estimates of the relevant statistical 
parameters' may increase the program's execution time unduly. 

(iii) One may abandon the statistical assessment of the alignments 
produced. This gives the greatest scope to input scoring systems, but 
precludes any reasonable scheme for automatic iteration of PSI- 
BLAST. Much experimentation will be required to detamine which 
of these approaches is the most finitfuL 

Realignment 

After the initial BLAST run, or a later PSI-BLAST iteration, the 
multiple aligimient used for subsequent iterations can be con- 
structed in a more sophisticated manner than described above. 
Rather than coalescing all the pairwise alignments that pass the 
threshold immediately into a multiple alignment, the most 
significant among them can be used to build an initial multiple 
alignment and associated position-specific score matrix, which 
can then be used to rescore and realign database sequences that 
received lower scores. This step can be iterated several times 
before another full-scale database search is executed. There are 
several potential advantages to this procedure, (i) Weaker 
pairwise alignments, that may be somewhat inaccurate, can be 
improved and perhaps extended before they are incorporated into 
the evolving multiple alignment, (ii) Unrelated sequences that 
received chance high scores can have their scores downgraded by 
an improved matrix, and perhaps be rejected before they are 
included in the alignment (iii) Related sequences that received 
relatively high alignment scores, but missed the threshold for 
inclusion, can have their scores increased, and perhaps be 
included in the multiple alignment. In short, the realignment 
procedure can prevent inaccurate pairwise aligrunents from 
corrupting the evolving multiple alignment, and can accelerate 
the recognition of related sequences, all for very little computa- 
tional cost. Preliminaiy studies suggest this line of developm^t 
to be a fruitful one. 

In conclusion, the new gapped version of BLAST is both 
considerably faster than the original one, and able to produce 
gapped alignments. While the relevant statistical parameters can 
no longer be calculated from theory, random simulation allows 
ftem to be estimated beforehand for commonly used amino acid 
substitution matrices and gap costs. For many queries, the 
PSI-BLAST extension can greatly increase sensitivity to weak 
but biologically relevant sequence relationships. PSI-BLAST 



retains the ability to report accurate statistics, per iteration runs in 
times not much greater than gapped BLAST, and can be used both - 
iteratively and fully automatically. These developments should 
enhance significantly thie utility of database search methods to the 
molecular biologist 

Note 

Source code for the new BLAST programs is available by 
anonymous ftp from the machine ncbijilmjiih.gov, within the 
directory 'blast', and the programs may be run from NCBFs web 
site at http://www.ncbi.nlm.nih.gov/ 
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